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SUMMARY 


A finite-difference solution code for the two-dimensional Navier-Stokes equa- 
tions has been combined with a moving-grid system. The thin layer Navier-Stokes 
equations with a turbulence model are solved in a time-accurate manner in order to 
study the unsteady aerodynamics around airfoils undergoing small amplitude pitching 
or heaving motions in the transonic regime. The accuracy of the solutions obtained 
by the use of the present moving-grid technique is investigated. The effects of the 
minimum grid size and the integrating time-step size on the solutions are also 
checked. Some of the solutions obtained by the present method are compared with 
experimental results. 

It is demonstrated that the unsteady aerodynamics around oscillating airfoils 
can be predicted fairly well by the present code for cases in which the dynamic 
angle of attack or displacement is small. 

INTRODUCTION 

Several schemes for the solution of unsteady gas dynamic equations have been 
developed in the last decade and some of these have been used to predict the 
unsteady aerodynamics around oscillating airfoils. 

As is well known, the unsteady flow fields can be simulated by using various 
levels of gas dynamic equations. At present, however, most of the numerical studies 
on unsteady aerodynamics are based upon the potential equation or the Euler equa- 
tions, that is, inviscid flow-field models. For example, Euler solutions have been 
obtained by Magnus and Yoshihara (1975) and Jameson (1975). The full-potential 
equation has been solved and studied by Isogai (1977, 1982, 1984), and Bridgeman and 
Steger (1982) in 2-D and 3-D problems. The small-disturbance equation has been 
studied by Goorjian and Guruswamy (1984). In their studies, the predictions of the 
unsteady aerodynamics agree well with experiments, especially when boundary layer 
corrections are employed. 

The inviscid unsteady transonic flow fields about oscillating airfoils can be 
determined at a low cost compared with obtaining the same information for viscous 
flow fields. So the solutions of the potential equation have already been applied 
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to aeroelasticity problems (Isogai, 1984; Bridgeman and Steger, 1982). There are a 
few studies of unsteady Navier-Stokes solutions: Chyu and Davis (1984) and Chyu and 

Kuwahara (1982) obtained Navier-Stokes solutions around oscillating airfoils. In 
their studies, the shock-induced separation was numerically simulated. Considering 
the capability of the recent super computers, however, it is time to apply the 
unsteady Navier-Stokes codes to practical problems. The viscous effects on tran- 
sonic airfoil response and stability are one of the major interests in aeroelas- 
ticity problems. 

To describe the effects of viscosity on oscillating airfoils, the Navier-Stokes 
equations must be solved correctly with an appropriate moving-grid system. The 
purpose of this report is to show the capability of the present code (based upon 
ARC2D) for solving unsteady problems, especially for oscillating airfoils in the 
transonic regime. 

In the present study, the unsteady aerodynamics around airfoils are computed by 
use of the thin-layer approximation. The mean angle of attack is zero and the 
amplitudes of the oscillations are small. In addition, the turbulence model is 

introduced in the computations with the assumption that the transition point is at 

the leading edge of the airfoil. 

Some of the results obtained by the present program are compared with those 
obtained by experiments. Although the amplitudes of the motions are small, it has 
been verified that the computed results are in good agreement with experimental 
results. The reduced frequencies of the cases studied here are relatively low, 
ranging from 0.1 to 0.4, which is based upon the full chord length. The computed 
airfoils are NACA 0012 and NACA 64A010. 

The author expresses his thanks to T. J. Barth for his offering numerical data 
about the hyperbolic grid around a NACA 0012 airfoil and J. Baeder for reading this 
manuscript. 

The author also expresses his thanks to the following Ames people: Dr. J. W. 

McCroskey for his kind advice and helpful discussions; Dr. T. H. Pulliam for his 

valuable comments; and Mr. M. Inouye for his kindness during the course of this 
study . 


THE OUTLINE OF SOLUTIONS 


The Navier-Stokes equations are written in a strong conservation-law form using 
Cartesian coordinates. We introduce the nondimensional quantities to scale the 
variables p (density), u,v (velocities), and e (the total energy) as 
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where ® indicates the quantities concerning freestream. Choosing a reference 
length, c, which is the airfoil chord length, time t may be nondimensional ized 
as t = ta^/c. The Reynolds number is expressed as 
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Note that Re uses a m and, therefore, Re based on u m (the usual case for 
experimentally given Reynolds number) must be scaled by M = u /a . For the 
remainder of this development, the ~ will be dropped for simplicity. 


Then we can write the Navier-Stokes equations in a nondimensionalized form as 


3, Q + 3 E + 3 F = Re ( 3 E + 3 F ) 
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Pressure is related to the conservative flow variables, Q, by the equation of 

state 
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(7) 


p = (y - 1 ) ^ p(u 2 + v 2 ) 

where y is the ratio of specific heats, generally set to be 1.4 for air. The 
speed of sound is given by a = yp/p for ideal fluids. The dynamic viscosity is 
m and is typically made up of a constant plus a computed turbulent eddy viscosity. 

The Navier-Stokes equations (Magnus and Yoshihara, 1975) may be transformed 
from Cartesian coordinates to general curvilinear coordinates assuming the relations 
between two coordinate systems as 


t = t 


5 = £(x,y,t) (8) 

n = n(x,y,t) 


We can easily obtain the relation between the Cartesian derivatives and the 
curvilinear derivatives by use of chain rule expressions as follows: 
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J is defined to be the metric Jacobian, that is, J - = (x^.y^ ~ ^y^)- Using 
the above equation and expressions, the Navier-Stokes equations are rewritten in a 
generalized curvilinear coordinate system 
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the contravariant velocities. The viscous flux terms are E □ J E + e F ) 

— 1 , » V XV V V 

and F =J ( n E + n F ) . 
v xv y v 


The stress terms, such as t xx are also transformed in terms of the E and n 
derivatives where 
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with terms such as u x expanded by chain rule. 

In high Reynolds number flows the effects of viscosity are concentrated near 
rigid boundaries and in wake regions. Most of the computing efforts should be 
concentrated on the resolution of the boundary layer. A grid is generated which is 
very fine in the normal direction to the body, and is coarse in the tangential 
one. Even if the full Navier-Stokes equations are solved on that kind of grid, the 
viscous terms associated with derivatives along the body will not be resolved. 


The thin-layer approximation requires that: 

1. All body surfaces be mapped onto coordinate surfaces. Specifically, 
n = constant coordinate surfaces. 

2. Grid spacing is clustered to the body surfaces such that sufficient resolu- 
tion for a particular Reynolds number is obtained. 

3. All the viscous derivatives in the E direction are neglected, while the 
terms in the n direction are retained. All of the inviscid terms are used. 
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Applying the thin-layer approximation to equations (12) and (15), the thin- 
layer Navier-Stokes equations are written 
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Applying an implicit three point time differencing scheme of the form, 
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(19) 


where AQ n = Q n+1 - Q n and Q n = Q(nAt). The parameters 9 and <)> can be chosen to 
produce different schemes of either first- or second-order accuracy in time. 

For 0 = 1 and $ = 0, we have the first-order Euler implicit scheme, and 
for 9 = 1 and <j> = 1/2, the three point implicit scheme, respectively. By solving 
equation (19) using the ADI (Beam-Warming) scheme, we can get the solution in a 
time-accurate manner. For more detail, readers are referred to Pulliam (1984). 


BOUNDARY CONDITIONS FOR UNSTEADY PROBLEMS 


We need boundary conditions for each conservation-law equation. At the solid 
surface, the density is extrapolated from the interior boundary to the flow field. 
For the momentum equations, the contravar iant velocity U and V are set to be zero 
at the body surface. At the trailing edge, the Kutta condition is applied for the 
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airfoil. The pressure is obtained by the normal momentum equation and the accelera- 
tion of the solid surface is taken into account as follows: 

' pU( V E + n ,V * ‘Vx * 5y"y ,P 5 * <n ? * °y )p n (20) 

When a C-type grid is used, a part of the boundary is split. All of the 
quantities on the cut line which are necessary for the next time step are evaluated 
by averaging the quantities just above and below the cut line. 

In the present study, the grid is changing during the unsteady computations: 
the inner boundary is moving as prescribed by a function of time although the shape 
of the airfoil is not changing (rigid-body motion). The far-field boundary is fixed 
with respect to the inertial coordinate system. The grid distribution between the 
inner boundary and the far-field boundary is changing. 

In the present technique, the orthogonality of the grid is not maintained even 
though the initial grid system is orthogonal. The orthogonality condition between 
£ = constant lines and the airfoil surface is important for the boundary conditions 
to be satisfied exactly. In our numerical experiments, however, in both cases of 
pitching motion and heaving motion, the violation of the grid orthogonality will not 
cause large errors or instability during the unsteady computations as long as the 
amplitude is small. 

At every time step, the velocity of each grid point is computed numerically 
using the coordinates at the previous time step and those of the updated grid. The 
accelerations of the grid points on the airfoil surface are also numerically calcu- 
lated in a similar manner. In the case of forced oscillations, the motion of the 
airfoil is prescribed as a function of time, so the velocities and accelerations can 
be computed analytically. However, in the present study, these quantities have been 
computed numerically for the purpose of checking the applicability of the program 
for the case of the numerical simulation of the aeroelasticity problems. In that 
case, the airfoil motions will not be prescribed in advance. As for the far-field 
boundary, uniform flow conditions are applied during the unsteady computations. 


MOVING-GRID SYSTEM 


For the solution of the unsteady Navier-Stokes equations, the grid should be 
changing with the airfoil motions. The boundary conditions must be satisfied on the 
moving surface at each time step. In the present study, a technique is introduced 
to update the grid system around an airfoil which is changing in location and atti- 
tude with respect to time. 

There are two techniques for a time varying-grid system: one is based upon the 

rigid motion of the grid system in the inertial coordinate system, and the other is 
based upon the deformation of the grid between the inner boundary and the outer 
boundary . 
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In the former technique, it is not necessary to update the instantaneous grid 
coordinates, and we need just metrics and their derivatives concerning space and 
time at each time step as long as we treat rigid airfoil sections. As to the 
computing time and storage, there are some advantages in this rigid motion of a grid 
system. Steger (1978), and more recently, Jesperson (1985), for example, solved 
Euler equations around a heaving airfoil in a time-accurate manner using this kind 
of grid. This technique, however, cannot be applied to problems in which the 
elastic deflections of airfoils should be taken into account. In 2-D problems, 
rigid-motion grid seems to be enough because the airfoil sections are considered to 
be rigid in most cases. 

On the other hand, in the latter technique, the inner boundary is changing with 
the motion of the airfoil while the outer boundary is at rest with respect to the 
inertial coordinate system. It will be applicable for aeroelasticity problems as 
well, although both the grid coordinates and the corresponding metrics must be 
updated at each time step. 

An implementation of the latter technique was presented by Chyu and Kuwahara 
(1982) and Chyu and Davis (1984). They successfully obtained the moving grid for 
the computations of pitching airfoils. First, they generated the static grids for 
the airfoil at the extreme angle of attack positions. The grid at each time step 
was obtained from those at the extreme positions by spatial interpolations. The 
inner boundary is changing with the motion of the airfoil and the far-field boundary 
is fixed in the inertial coordinate system. In their technique, at least two grid 
systems must be stored, even in the simple case of a pitching airfoil. 

In the present study, a method is proposed to generate a time-varying grid 
system. In this method, the grid points near the inner boundary are moving com- 
pletely with the motion of the airfoil while the far-field boundary is at rest 
during the unsteady computations. The grid distribution between the airfoil surface 
and the far-field boundary is determined by an analytical function along the 
5 = constant lines. From the viewpoint of the accuracy of the computations, the 
deformations of the grid should be minimum in the flow field. The analytical func- 
tion is determined so that the deformations of the grid in the vicinity of the 
airfoil are small, but those in the neighborhood of far-field boundary are rela- 
tively large. In problems with elastic deflections, some special technique must be 
employed to specify the grid points on the airfoil surface. Once the grid distribu- 
tions on the rigid or elastic airfoil are specified, the grid distributions in the 
field will be determined by use of an algebraic function as 


Ax.^ = Ax. • Vi - (n*) 2 

(21) 

Ay i,j = Ay i ‘ ^ " (n * )2 

(22) 
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where 


n* = 0 ; n < n Q 


ti* = (n - n )/( n - n ) ; n > n 
o max o ’ - o 


The increments of the grid coordinates are denoted by Ax^ • and Ay^ ■ , and 
Ax i and Ay i are the displacements of the grid points on the airfoil (and the cut 
line when a C-type of grid is used) in the x and y directions, respectively. A 
parameter n 0 takes a value between 0 and n max depending upon the airfoil motion 
of interest. Usually, a value between 0.6 and 0.7 is suitable for n Q . The vector 
(Axj, Ay^ is computed from the analytical function of t in cases of forced oscil- 
lations or obtained by use of the equations of motion of the airfoil in cases of 
aeroelasticity problems. The updated grid will be 
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The grid velocities and accelerations are computed by use of Ax^ j and Ay^ j 
and the new coordinates can overwrite old ones. 

Other choices for updating the grid system are possible. The grid has to be 
changed at each time step. The computing time to update the grid system should be 
small, and it is desirable that the scheme of changing the grid be simple and reli- 
able for any kind of airfoil motions. 

To demonstrate the present technique, a moving grid-system is shown. The 
airfoil is rigid and its motion around the quarter-chord is described as 

h = 0.0 + 0.1 • sin(wt + 30°) (25) 

a = 0.0 + 10° • sin(wt) (26) 

that is, a heaving and pitching coupled motion with amplitudes, 0.1 and 10.0 deg, 
respectively. Both of the amplitudes are large in order to illustrate the useful- 
ness of the present method. The phase difference between heaving and pitching 
motions is 30.0 deg. The initial grid is a hyperbolic-C-type generated by use of 
the method of Kinsey (1984). 

Unfortunately, the orthogonality condition between the body surface and the 
^-constant lines is violated in spite of the fact that the initial grid is generated 
so as to satisfy orthogonality. However, the grid at each time step looks very 

reasonable. At every ut = 30°, the grid is shown in figure 1. The grid motion is 
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very smooth even though both of the amplitudes of the heaving and pitching motions 
are relatively large. 

The C-type of grid has a cut line extending from the trailing edge to the 
downstream boundary. At each time step, the cut line is reset smoothly so as to 
approximate a streamline. After setting the new cut line, the grid distribution in 
the region behind the trailing edge of the airfoil is specified in the same manner 
as prescribed above. 

The computing time to update the grid at each time step is less than 1 % of that 
for solving the flow field, including the recalculation of grid coordinates, grid 
velocities, grid accelerations, and other quantities necessary to compute the flow 
field. In the present technique, the grid pattern is deforming with the flow 
field. Thus, the Jacobians and the metrics at the grid points are also changing 
with respect to time. The changing Jacobians might have some effect on the flow 
computations since the present code does not take them into account. We, therefore, 
have to check the effect by comparing two solutions: one obtained by using a 

deformed moving grid, the other by using a rigid body motion grid. The results of 
our comparison will be shown below. 


RESULTS 


Some verifications of the computed results are presented. Throughout the 
present study, the computations are performed under the flow condition 
Re = 12x10 . The thin-layer approximation and the Baldwin-Lomax turbulence model 
are used assuming the transition point is at the leading edge of the airfoil. 
Uniform flow conditions are applied at the far field boundaries which are located 
ten chord lengths from the airfoil. The grid size is 159 x 51. The grids used in 
these studies are elliptic-C-type , except for those used to demonstrate the moving 
grid technique. Thus, the grid for the steady state (that is, the initial grid for 
unsteady computations) is generated by solving a set of nonlinear Poisson's equa- 
tions. The minimum size of grid and the way of stretching the grid from inner 
boundary to outer boundary is chosen by specifying parameters in the equations of 
Sorensen (1980). It satisfies the orthogonality condition between the airfoil 
surface and the ^-constant lines, though some skewness is observed, especially 
along the cut line after the trailing edge. 

The steady state is first computed. After the steady solution is sufficiently 
converged, the airfoil is forced to oscillate periodically. For the steady solu- 
tions, the Euler-Implicit scheme is employed; while for unsteady solutions, the 
Three-Point-Implicit scheme is used. The airfoil motion and the flow field around 
it are simulated for three cycles. In the third cycle, the data of interest are 
stored. The quantities obtained in a time history are analyzed by Fourier analy- 
sis. The definitions of the steady state, in-phase component, and out-of-phase 
component are 
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respectively. In the above expressions, F(t), is an arbitrary function of time 
and 6 is the amplitude of the airfoil oscillation, and k is the reduced fre- 
quency defined by 



co 


(30) 


The major capabilities of the present program to be investigated are: 

1. The accuracy of the computations using the moving grid system. 

2. The effects of the minimum grid size used for the unsteady computations. 

3. The effect of the time integration step size. 

4. Comparisons between the computational results and experimental data. 

It is actually difficult to watch the pure effects of each of the above 
items. For example, the magnitude of the smoothing parameters which is necessary to 
maintain the stability of the calculations will depend on the time integration step 
size and the minimum grid size. However, it is useful to check the above items for 
the numerical studies on the unsteady viscous aerodynamics. 

Fortunately, for the case of an airfoil oscillating in a rigid motion, we can 
alternatively use the rigid motion grid method. The Jacobian at each grid point is 
changing with respect to time when the present moving-grid technique is employed. 
However, in the case of the rigid motion grid method, it is constant. The accuracy 
of the solutions obtained using the present technique can be checked by comparing 
the two results. Two computed results are compared to check the effects of deforma- 
tion of the grid using a NACA 0012 airfoil in a pitching motion around the midchord 
axis. The amplitude of the motion is 1.0° and the reduced frequency is 0.2 with a 
freestream Mach number of 0.8. 
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As the C time histories of both cases are compared, no differences can be 
recognized. The in-phase and the out-of-phase components for both cases are com- 
pared in figure 2(a), (b). In this figure, the results from the present method are 
illustrated by solid lines and those from the rigid motion grid are illustrated by 
dotted lines. There are no significant differences between both results in the 
calculations of in-phase and out-of-phase components in front of the shock wave; 
but, at the peak of the shock wave and behind the shock wave, there are slight dis- 
crepancies between the two results. 

In the computation using the rigid motion grid, some special care was taken 
into account concerning the evaluation of the grid velocities. Fortunately, in the 
case of pitching motion, the grid velocity is computed as accurately as possible 
using analytical relations between the inertial coordinate system (x,y) and the 
coordinate system fixed on the airfoil (x*,y*). We can easily obtain the following 
relations when the pitching axis is chosen as the origins of both systems. 
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r 2 = x 2 + y 2 = (x*) 2 + (y*) 2 and 0 
6 ( t ) = 6 sin t(w • t) 
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the airfoil given by 
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(34) 


By use of the above relations, the calculation can be continued during the pitching 
motion without recalculating the grid coordinates. 

In figure 3(a), (b), the change in the solutions caused by the minimum grid 
sizes are illustrated; one is by use of minimum grid size, 2.0x10” and the other is 
by 2.0x10”^. The former is shown by solid lines and the latter shown by dotted 
lines. The experimental results are also plotted in the figures by triangles. The 
airfoil computed here is NACA 64A010 in a pitching motion around a quarter-chord 
axis with reduced frequency of 0.4 and amplitude of 1.0°. The freestream Mach 
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number is 0.8. The results are similar with slight differences observed ahead of 
the shock wave. 

According to our experiences, it is sometimes difficult to get the steady-state 
solutions sufficiently converged when a fine grid is used. It is especially true 
with the computations of the NACA 64A010 airfoil using an elliptic-C-type grid. The 
unsteady aerodynamics will depend upon the convergence status of the steady-state 
solution which becomes the initial condition for the unsteady problem. Thus, the 
steady-state solution must be examined before starting unsteady computations. 

The minimum grid size has a close correlation with the number of grid points in 
the boundary layer. The differences between the two solutions are due to the reso- 
lution of the boundary layer. It is clear from the present results that the finer 
the grid used, the better is the agreement of computational solutions with experi- 
mental data. But, it should be noted that the accuracy of the solutions will 
depend, not only upon the minimum grid size, but also upon the manner of stretching 
the grid, especially in unsteady problems. The minimum grid size of 2.0x10“° is 
sufficient to obtain good agreement with experiments. In all of the remaining 
results, the minimum grid size of 2.0x10“° is used. 

The effect of the integrating time step size is also investigated. The two 
results are shown in figure 4(a), (b). The airfoil is a NACA 0012 which is in a 
pitching motion with reduced frequency, 0.1628. The mean angle of attack is 0° and 
the amplitude of the pitching motion is 2.51°. The Mach number is 0.755. One is 
the result obtained by use of the time integrating step size, 0.01692; the other is 
by use of a time integrating step size of half that size, that is, 0.00846. They 
are based upon the nondimensionalization of t = ta m /c and each time step size 
results in 3020 and 6040 iterations per cycle, respectively. The former is shown by 
the solid line and the latter is shown by the circles. The solutions are almost 
identical. For the unsteady problems, it is inevitable to adjust the smoothing 
parameters in the calculations because of stability reasons. In the present study, 
we used the minimum value of the smoothing parameters to maintain stability of the 
unsteady computations. The values of the smoothing parameters should be as small as 
possible. However, for practical computations, there will be a trade-off between 
the accuracy and the stability. 

Figure 5 demonstrates the unsteady aerodynamics of a NACA 64A010 airfoil in a 
heaving motion. The reduced frequency is 0.4 and freestream Mach number is 0.8. 

The amplitude of the heaving motion is 0.0436 based upon the full chord length. The 
same problem was computed by Steger (1978) and by Jesperson (1985) using an Euler 
code. 


In figure 6(a), (b), unsteady aerodynamic components are compared with experi- 
mental results. The case is a NACA 64A010 airfoil in a pitching motion. The 
reduced frequency is 0.1 and the pitching axis is a quarter chord point. The ampli- 
tude of the motion is 1.0° with a freestream Mach number, 0.8. In the figures, the 
experimental results are shown by circles and the computed results by solid lines. 

In figure 3, the solutions for the case of reduced frequency, 0.4, were already 
compared. In most of the experimental studies, the number of pressure orifices is 
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insufficient and the peak of the unsteady pressure at the shock wave cannot be 
captured precisely. Thus, it is somewhat difficult to compare the unsteady compo- 
nents of the pressure near the shock wave. As long as this fact is taken into 
account, both components agree fairly well with the experimental data, except for 
the out-of-phase component in the shock wave region. The computed results show a 
peak at the shock wave in the out-of-phase component while the experiment data show 
just a step change there. In figure 3, the same trend is observed. That is the 
major difference between the computed results and experiment. 

Finally, a NACA 0012 airfoil is studied. The case is a NACA 0012 airfoil in a 
pitching motion around a quarter-chord axis. The freestream Mach number is 0.755 
and the reduced frequency is 0.1628. The amplitude of the pitching motion is 2.51°, 
which is relatively large. In figure 7, the steady state solution is shown. The 
shock wave is very weak and smeared. This result compares well with the results 
presented in Goorjian and Guruswamy (1984). Starting with this steady-state solu- 
tion, three cycles are simulated as was done in the other problems. Figure 8 shows 
a comparison at every 45° during the third cycle of the pitching motion. The com- 
puted results are shown by solid lines and the experimental data are plotted by 
daggers (t). The mean angle of attack and the Reynolds number are 0.0 and 12x10^ in 
the computation; they are 0.016 and 5.5x10^ in the experimental data (Lambourne, 
1985). For this large amplitude of the pitching motion, the shock wave is rela- 
tively strong in some range of phase angle, and the unsteady component before the 
shock wave is very large. Thus, both components of the unsteady C distribution 
indicate no peaks at the shock wave, but present just step change at the shock 
location. In figure 8, most of the computed results are in excellent agreement with 
the experimental results, except for the part just after the shock wave. This is 
due to the lack of the capability of resolving the shock and boundary layer 
interactions. 

The time history of the Mach number contours around the airfoil is illustrated 
in figure 9. The number of grid points is insufficient, especially around the shock 
wave. The interaction between the shock wave and the boundary layer cannot be 
observed precisely. Qualitatively, however, we can see the change of the flow field 
around airfoil in a cycle of motion. In figure 10, the unsteady part of the aerody- 
namics was obtained by Fourier analysis. However, the number of experimental data 
points is too small for Fourier analysis. Thus, the experimental results are not 
presented in the figure. 

As for the sensitivity of solutions to grids, it is difficult to reach a con- 
clusion because no single parameter can be isolated. If the grid is changed, there 
might be some effects even on the steady solutions. In the present study, just some 
of the tendencies of the solutions are shown concerning the minimum grid size as a 
result of number of numerical studies on unsteady Navier-Stokes solutions. 
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CONCLUSIONS AND DISCUSSIONS 


A moving-grid technique has been demonstrated. Some of the obtained solutions 
have been compared with experimental results. 

The present moving-grid technique is very powerful in the case of small ampli- 
tudes, and it does not increase the total computing time significantly. In the 
present study, both hyperbolic-C-type and elliptic-C-type grids were used. However, 
this moving-grid technique is applicable to other kind of grids. It has been shown 
that the effect of the deformations of the grid is very small and the accuracy of 
the solution can be retained during the unsteady computations. 

The effect of minimum grid size has also been investigated. The best results 
were obtained when the minimum grid size of 2x10 - ^ was used. Cases may differ, but 
generally speaking, the finer the grid used, the better the results that are 
obtained. This is because the boundary layer should be resolved as accurately as 
possible, especially in the unsteady problems. It cannot be emphasized too strongly 
that most of the computing efforts should be concentrated on the resolution of the 
boundary layer. Some of the solutions have been shown in this report. However, it 
must be stressed that even the solutions which look reasonable and compare well with 
experimental results need more investigation concerning the boundary-layer behavior, 
especially for problems with shock waves. 

The effect of time-integration step size has also been investigated. To obtain 
unsteady solutions of the Navier-Stokes equations in a finite and realistic comput- 
ing time, some smoothing terms are inevitable for stability reasons. It is defi- 
nitely better to use as small smoothing parameters as possible, even though a very 
small time step must be used for time integrations. As mentioned previously, it is 
a trade-off between accuracy and computing time. Fortunately, according to the 
present study, the solutions are independent of the time-integration step size as 
long as they are small enough. But, it is noted that when the time-integration step 
size is changed, the smoothing parameter is also adjusted to maiantain stability or 
accuracy. 

In the Navier-Stokes computations, there are many parameters to be determined 
before computations. The mathematical model for the computations is determined when 
all of the parameters are fixed, including the grid system. The manner of stretch- 
ing the grid has some effects on the unsteady solutions. In the viscous flow field, 
there are large dissipations. However, even an inner grid line can be a boundary 
for some disturbances. If the stretching is poorly done, certain waves cannot 
propagate through the inner grid lines and might be reflected as if from a rigid 
boundary. In that case, the solution will never become periodical, even if the 
motion of the airfoil is periodic. 

The computed results were compared with experimental ones in some cases. The 
obtained results are more than satisfactory, but there are some disagreements about 
the shock wave motions. It is due to the lack of grid points around the shock wave, 
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and the fact that the shock wave-boundary layer interactions are not resolved ade- 
quately. Moving adaptive grid should be introduced for further investigations. 

The conclusions of the present study are: 

1. The accuracy of the computations can be maintained during unsteady computa- 
tions using the present moving-grid technique. 

2. Although it depends on requirements of accuracy, the minimum-grid size 
should be of the order of 10”^. (Of course, it will depend on the Reynolds number.) 

3. The unsteady solutions are independent of the time-integration step size as 
long as it is small enough. 

4. The computed results are reasonable and compare well with experimental data 
in cases of small amplitude. 

For further investigations, the following items should be addressed: 

1. The sensitivity of the unsteady solutions to grids; the effects of distance 
to the far-field boundary and the grid stretching pattern. 

2. The frequency response of the aerodynamics produced by the present program 
(especially in the case of small amplitude). 

3. The prediction capability of shock and boundary layer interactions. 

4. Computations of cases of large amplitudes. 

Item 2 directly above is closely related to item 1. For practical applica- 
tions, we need to develop a code which includes a systematic grid generator. Either 
of items 3 or 4 may be another major problem beyond the present study. 
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Figure 2.- Comparisons between present grid and rigid motion grid, (a) In-phase 
(NACA 0012, = 0.8, 6a = 1.0°, k = 0.2). (b) Out-of-phase (NACA 0012, 

M m = 0.8, 6a = 1.0°, k = 0.2). 



Figure 3.- Effects of minimum grid size on solutions, (a) In-phase (NACA 64A010, 
= 0.8, 6a = 1.0°, k = 0.4). (b) Out-of-phase (NACA 64A010, M^, = 0.8, 

6a = 1.0°, k = 0.4). 


21 






Figure 4.- Effects of time integrating step size on solutions, (a) In-phase (NACA 
0012, M ro = 0.755, «a = 2.51°, k = 0.1628). (b) Out-of-phase (NACA 0012, 

M, = 0.755, 6a = 2.51°, k = 0.1628). 



Figure 5.- Unsteady aerodynamics: NACA 64A010 in heaving (M m = 0.8, 6h = 0.0436, 

k = 0.4). 
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Figure 6.- Unsteady aerodynamics: NACA 64A010 in pitching, (a) In-phase (M^ = 0.8, 

X a = 0.25, So □ 1.0°, k = 0.4). (b) Out-of-phase (M m = 0.8, X a = 0.25, 

6a = 1.0°, k = 0.4). 


NACA0012 STEADY STATE 

159X51 GRID, M = 0.755, a =0.0°, Re=1.20X10 7 

' OO ' ' 



Figure 7.- Steady state solution: NACA 0012 (M m = 0.755, a m = 0.0°) 




Figure 8.- Time history of C distributions: NACA 0012 (M m = 0.755, X a = 0.5, 

6a = 2.51°, k = 0.1628). 
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Figure 8,- Concluded. 
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MACH NUMBER CONTOURS 

NACA0012 IN PITCHING MOTION (K = 0.1628, XA = 0.25) 
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Figure 9.- Concluded. 
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Figure 10.- Unsteady aerodynamics: NACA 0012 in pitching (M m = 0.755, X 

Sa = 2.51°, k = G..1628). 


1. Report No. 2. Government Accession No. 

NASA TM-8834 1 

3. Recipient's Catalog No. 

4 Title and Subtitle 

A VERIFICATION OF UNSTEADY NAVIER-STOKES SOLUTIONS 
AROUND OSCILLATING AIRFOILS 

5. Report Date 

September 1986 

6. Performing Organization Code 

ATP 

7. Author(s) 

Jiro Nakamichi 

8 Performing Organization Report No. 

A-86337 

10. Work Unit No. 

9 Performing Organization Name and Address 

Ames Research Center 
Moffett Field, CA 94035 

1 1 . Contract or Grant No. 

13 Type of Report and Period Covered 

Technical Memorandum 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, DC 20546 

14 Sponsoring Agency Code 

505-60 

1 15 Supplementary Notes 


Point of contact: Mamoru Inouye, Ames Research Center, MS 202A-1, Moffett 

Field, CA 94035 (415) 694-5126 or FTS 464-5126 


16 Abstract 

A finite-difference solution code for the two-dimensional Navier-Stokes 
equations has been combined with a moving-grid system. The thin layer 
Navier-Stokes equations with a turbulence model are solved in a time- 
accurate manner in order to study the unsteady aerodynamics around airfoils 
undergoing small amplitude pitching or heaving motions in the transonic 
regime. The accuracy of the solutions obtained by the use of the present 
moving-grid technique is investigated. The effects of the minimum grid 
size and the integrating time-step size on the solutions are also checked. 
Some of the solutions obtained by the present method are compared with 
experimental results. 


It is demonstrated that the unsteady aerodynamics around oscillating 
airfoils can be predicted fairly well by the present code for cases in 
which the dynamic angle of attack or displacement is small. 


17. Key Words (Suggested by Author(s)) 


18. Distribution Statement 



Oscillating airfoils 
Thin layer Navier-Stokes 
Moving-grid technique 

solutions 

Unlimited 

Subject category - 02 

19 Security Classif. (of this report) 

20. Security Classif. (of this page) 

21. No. of Pages 

22. Price" 

Unclassified 

Unclassified 

30 

AO 2 


For sale by the National Technical Information Service, Springfield, Virginia 221 61 




















